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Recent developments in the numerical renormalization group (NRG) allow the construction of 
the full density matrix (FDM) of quantum impurity models (see A. Weichselbaum and J. von Delft 
in Ref. 1) by using the completeness of the eliminated states introduced by F. B. Anders and 
A. Schiller in Ref. 2. While these developments prove particularly useful in the calculation of 
transient response and finite temperature Green's functions of quantum impurity models, they may 
also be used to calculate thermodynamic properties. In this paper, we assess the FDM approach 
to thermodynamic properties by applying it to the Anderson impurity model. We compare the 
results for the susceptibility and specific heat to both the conventional approach within NRG and to 
exact Bethe ansatz results. We also point out a subtlety in the calculation of the susceptibility (in a 
uniform field) within the FDM approach. Finally, we show numerically that for the Anderson model, 
the susceptibilities in response to a local and a uniform magnetic field coincide in the wide-band 
limit, in accordance with the Clogston- Anderson compensation theorem. 

PACS numbers: 75.20.Hr, 71.27.+a, 72.15.Qm 



I. INTRODUCTION 

The numerical renormalization group method, 3-6 has 
proven very successful for the study of quantum impu- 
rity models. . Initially developed to describe, in a con- 
trolled non-perturbative fashion, the full crossover from 
weak to strong coupling behavior in the Kondo problem 3 
and the temperature dependence of the impurity ther- 
modynamics ; 3 ~ 5 > 8 it has subsequently been extended to 
dynamic 9-11 and transport properties 12 of quantum im- 
purity models. Recently, a number of refinements to the 
calculation of dynamic properties have been made, in- 
cluding the use of the correlation self-energy in evaluating 
Green functions 13 , the introduction of the reduced den- 
sity matrix , and the introduction of a complete basis 
set using the eliminated states in each NRG iteration 2 . 
The latter, in combination with the reduced density ma- 
trix, has been used to evaluate the multiple shell summa- 
tions arising in the time dependent transient response in 
quantum impurity problems 2 ' 15 and offers the possibility 
to investigate truly non-equilibrium steady state trans- 
port within the NRG method 16 . In addition, the com- 
plete basis set offers an elegant way to calculate finite 
temperature Green functions which satisfy the fermionic 
sum rules exactly 1 ' 17,18 For recent applications of this 
technique to transport properties, see Refs. 19 and 20. 

In this paper we benchmark the full density matrix 
(FDM) approach to thermodynamic properties, by ap- 
plying it to the prototype model of strong correlations, 
the Anderson impurity model 21 . This model has been 
solved exactly using the Bethe ansatz. 22-28 A numerical 
solution of the resulting thermodynamic Bethe ansatz 
(TBA) equations therefore allows one to compare the 
FDM results for quantities such as the specific heat and 



the susceptibility with essentially exact calculations from 
the Bethe ansatz. In addition, we shall also compare the 
FDM results for specific heats and susceptibilities with 
those of the conventional approach 29 (see Sec. II for a 
more precise definition of what we term "conventional" ) . 



The paper is organized as follows. In Section II we 
specify the Anderson impurity model, and outline the 
conventional approach to thermodynamics within the 
NRG method 29 . The FDM approach to thermodynam- 
ics is described in Section III. Section IV contains our 
results. The impurity contribution to the specific heat, 
Ci, np , calculated within the FDM approach, is compared 
with Bethe ansatz calculations 22-28 and to calculations 
using the conventional approach in Sec. IV A. The im- 
purity contribution to the susceptibility, Ximp, with the 
magnetic field acting on both impurity and conduction 
electron states, calculated within FDM is compared with 
corresponding results from Bethe ansatz and the con- 
ventional approach within NRG in Sec. IV B). Results 
for the Wilson ratio as a function of the local Coulomb 
repulsion and the local level position are also given in 
Sec. IV B). In Sec. IV C we consider also the local sus- 
ceptibility of the Anderson model, xioo with a magnetic 
field acting only on the impurity, and show by compari- 
son with Bethe ansatz results for Ximp, that xioc = Ximp 
for both the symmetric and asymmetric Anderson model 
in the wide-band limit. In addition, we also compare 
the FDM and conventional approaches for another local 
quantity, the double occupancy, in Sec. IV D. Section V 
contains our summary. Details of the numerical solution 
of the thermodynamic Bethe ansatz equations may be 
found in Ref. 30. 
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II. MODEL, METHOD AND CONVENTIONAL 
APPROACH TO THERMODYNAMICS 

We consider the Anderson impurity model 21 in a mag- 
netic field B, described by the Hamiltonian 

H = H imp + H o + H int + H B . 

The first term, H lmp = J2 a edd^da + Un d ^n di , de- 
scribes the impurity with local level energy e d and on- 
site Coulomb repulsion U, the second term, Hq = 
Sfecr e fc c lo- c fco-i is the kinetic energy of non-interacting 
conduction electrons with dispersion e^, the third term, 
H-mt = VY^ka( c kad<? + d\c ka ), is the hybridization be- 
tween the local level and the conduction electron states, 
with V being the hybridization matrix element, and, 
the last term Hb = —gfJ-BBS z .tot where S z jot is the 
z-component of the total spin (i.e. impurity plus conduc- 
tion electron spin), is the uniform magnetic field acting 
on impurity and conduction electrons, g is the electron 
g-factor, and fiB is the Bohr magneton. We choose units 
such that g = fj,^ = 1 and assume a constant conduction 
electron density of states per spin N(e) — 1/2D, where 
D = 1 is the half-bandwidth. The hybridization strength 
is denoted by Ao = ttV 2 N(0) and equals the half-width 
of the non-interacting resonant level. 

The NRG procedure 4-1 ' consists of iteratively diago- 
nalizing a discrete form of the above Hamiltonian H. 
It starts out by replacing the quasi-continuum of con- 
duction electron energies —D < < D by logarithmi- 
cally discretized ones about the Fermi level Ef = 0, i.e. 
e„ = ±DK~ n ~( x ~ z > , n — 1, . . . , where A > 1 is a rescal- 
ing factor. Averaging physical quantities over several re- 
alizations of the logarithmic grid, defined by the param- 
eter z € (0,1], eliminates artificial discretization induced 
oscillations at A ^ i 29,31,32 Rotating the discrete con- 
duction states into a Wannier basis f n(T ,n = 0, 1, 2, . . . 
at the impurity site, one arrives at the form H = 
lim m _ >00 H m , where the truncated Hamiltonians H rn are 
defined by H rn = H imp + H hyh + YZ =0cr £n( z )fno-fna + 
Er=ol i «( z )(/n ( r/n+i<T + fl+iafna), with H imp as de- 
fined previously and H hyh = VY,ML d <* + d tfoo-)- Tnc 
onsite energies e„(z) and hoppings t n {z) reflect the en- 
ergy dependence of the hybridization function and den- 
sity of states. 4-6 The sequence of truncated Hamilto- 
nians H m is then iteratively diagonalized by using the 
recursion relation H. m+1 = H m + J2 a ^m{^)Smafma + 

E<t tmWiflurfm+ia + fL+Xcrfm*)- Tne resulting eigen- 
states \p; to) and eigenvalues E™ , obtained on a decreas- 
ing set of energy scales uj m (z) ~ t m (z), m = 0, 1, . . . , are 
then used to obtain physical properties, such as Green's 
functions or thermodynamic properties. Unless other- 
wise stated, we use conservation of total electron number 
N e , total spin S, and total z-component of spin S z in the 
iterative diagonalization of H at B = 0, so the eigenstate 
\p;m) is an abbreviation for the eigenstate \N e SS z p;m) 
of H m , with energy E^^ Sp (abbreviated as E™) where 
the index p = 1, . . . distinguishes states with the same 



conserved quantum numbers. As long as m < too — 1, 
where typically too = 4 — 6 , all states are retained. 
For to > m , only the lowest energy states are used to 
set up the Hamiltonian H m+ \. These may be a fixed 
number iVk OC p of the lowest energy states, or one may 
specify a predefined mo, and retain only those states 
with rescaled energies (E™ — Eq S ) /t m (z) < e c (A), where 
Eq S is the (absolute) groundstate energy at iteration to 
and e c (A) is A-dependent cut-off energy. 31 ' 33,34 For most 
results in this paper, we used too = 4, 5, respectively 
for H, H , and e c (A) = 15-y/A « 47 for A = 10 and 
found excellent agreement with exact continuum results 
from Bethe ansatz (after appropriate z-averaging, see be- 
low). Calculations at smaller A = 4, using too = 5,6 
for H, H , respectively, and e c (A) = 40 were also car- 
ried out for the local susceptibility in Sec. IV C. These 
showed equally good agreement with corresponding con- 
tinuum Bethe ansatz results, indicating that the too and 
A-dependence of our results (after z-averaging) is negligi- 
ble. In our notation, the number of retained states at it- 
eration to (before truncation) grows as 4 m+1 , so the value 
of to cannot be increased much beyond 5, in practice, 
due to the exponential increase in storage and computer 
time. As our calculations show, this is also not necessary, 
since agreement with exact Bethe ansatz calculations is 
achieved already for to > to = 4, 5. 

The impurity contribution to the specific heat is de- 
fined by Cimp(T) = C(T)-C (T), where C(T) and C (T) 
are the specific heats of H and Hq, respectively. Simi- 
larly, the impurity contribution to the zero field suscep- 
tibility is given by x imp (T) = x(T) - Xo(T), where x( T ) 
and Xo(T) are the susceptibilities of H and Hq, respec- 
tively. Denoting by Z(T,B) and Z a (T,B) the partition 
functions of H and Hq, and fl(T, B) = -k B TlnZ(T, B) 
and n (T,B) = -k B T\nZ (T,B) the corresponding 
thermodynamic potentials, we have, 

C ( T ) = = ^((H < H)f), (1) 

Ca(T) = -T 92 °0± =k B f({HQ-(H )Y)Q, (2) 

x(T) = - d2 f£ 2 B) \ B =o = P{gHB?(Sl tot ), (3) 

Xo(T) = - dH f§ B ± \ B = = P{g^) 2 (Sl tot \, (4) 

where S ztot is the z-component of total spin for Hq. 

We follow the approach of Ref. 29, which we term 
the "conventional" approach, to calculate the thermo- 
dynamic averages appearing in (1-4) at large A ^> 1 
(thermodynamic calculations at smaller values of A < 3 
are also possible, 4 ' 12 however, truncation errors increase 
with decreasing A): for any temperature T, we choose 
the smallest to such that kgT > t m (z) and we use 
the eigenvalues of H m to evaluate the partition func- 
tion Z m (T) = J2 P e ~ E ™^ kBT an d the expectation val- 
ues appearing in (1-4). Calculations for several z = 
(2i — l)/2n z ,i = l,...,n z with typically n z = 2,4 or 
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8 are carried out and averaged in order to eliminate dis- 
cretization induced oscillations at large A 3> 1. A dense 
grid of temperatures defined on a logarithmic scale from 
10~ 4 Tk to 2D was used throughout, where Tk is the 
Kondo scale for the symmetric Anderson model for given 
U [see Eq. (20)]. An advantage of the FDM approach, 
which we describe next, is that such a dense grid of tem- 
peratures can be used without the requirement to choose 
a best shell for a given T and z. This is possible within 
the FDM approach, because the partition function of the 
latter contains all excitations from all shells. 



III. THERMODYNAMICS WITHIN THE FDM 
APPROACH 

An alternative approach to thermodynamics is offered 
by making use of the eliminated states 2 from each NRG 
iteration. These consist of the set of states \lem) — 
|Zm)|e) obtained from the eliminated eigenstates, \lm), 
of H m , and the degrees of freedom, denoted collectively 
by e, of the sites i = m + 1, N, where N is the 
longest chain diagonalized. The set of states \lem) for 
m = m , . . . , N form a complete set, with completeness 
being expressed by 2 



N 



1= J2\lem')(lem'\ 



(5) 



m'—mo 



where mo — 1 is the last iteration for which all states 
are retained. Weichselbaum and von Delft 1 introduced 
the full density matrix (FDM) of the system made up of 
the complete set of eliminated states from all iterations 
to = toq + 1, . . . , N. Specifically, the FDM is defined by 



P : 



E £i< 

m=m le 



-PET 



em) 



Z{T) 



(lem\ 



(O 



where Z(T) is the partition function made up from the 
complete spectrum, i.e. it contains all eliminated states 
from all H m , to = too, . . . ,N [where all states of the last 
iteration m = N are included as eliminated states, so 
that Eq. (5) holds]. Similarly, one may define the full 
density matrix, po, for the host system Hq, by 



JV 



Po 



E Ei' 



-PET 



em) 



(lem\ 



(7) 



where Zq(T) is the full partition function of Hq. Note 
that m may differ from too, as the impurity site is miss- 
ing from Hq- In order to evaluate the thermodynamic 
average of an operator O with respect to the FDM of 
equation (6), we follow Weichselbaum and von Delft 1 
and introduce the normalized density matrix for the m'th 
shell in the Hilbert space of Hm'- 



-PEV 



(le 



(8) 



Normalization, Tr[/5 m ] = 1, implies that 



i = E 



,-PEV 



jN—m Zm 

a . , 

^m. 



(9) 



where Z m = J2i e~P E ™ and Z m = d N - ,n Z m with the 
factor d N ~ m resulting from the trace over the N — m en- 
vironment degrees of freedom e = (e m+ i, e m+ 2, . . . , e^r). 
For the single channel Anderson model, considered here, 
d = 4, since each assumes four possible values (empty, 
singly occupied up/down and doubly occupied states). 
Then the FDM can be written as a sum of weighted den- 
sity matrices for shells to = mo , . . . , N 



N 

p= E Wm p~ m ' 

m—mo 

y 

m -a z , 



(10) 



(11) 



where „ Wm = 1 and the calculation of the weights 
w m is outlined in Ref. 20. 

Substituting p = Y^m' Wm ' Pm' m ^ the expression for 
the thermodynamic average (O) and making use of the 
decomposition of unity Eq. (5), we have 



(6) p = Tr [pO 



-PET 



(l'e'm'\0 w m \lem) — -= (lem\l'e'm'} 



V e' m' 



J2o 



lem 



II w m - 



(12) 



Y,d N - m w m O m 

lm 
N 

E W rnO\f— 

1 > 

m-mo,i 



e -pEr 
n ^ u d N - m z„ 



(13) 



where orthonormality (lem\l'e'm') = 5iir6 ee /5 mm i, and 
the trace over the N — m environment degrees of free- 



dom J2i. 



= Elm d 



N- 



. . has been used and = 
(lm\0\lm) . A similar expression applies for expectation 
values (O) po with respect to the host system Ho. For each 
temperature T and shell m, we require w m (T) and the 
factor BJ"-(T) = e~ pE ™ /Z m where Z m = J2i e~ PE >" ■ Nu- 
merical problems due to large exponentials are avoided 



by calculating B" 1 (T) = e 



-p(ET 



>o )/Z' m where Z' m 



e@ E ™ Z m and E™ is the lowest energy for shell m. 

The calculation of the full partition function Z(T) , like 
the weights w m (T), requires care in order to avoid large 
exponentials (see Ref. 20). Note also, that the energies 
EJ"" in the above expressions denote absolute energies of 
H m . In practice, in NRG calculations one defines rescaled 
Hamiltonians H m in place of H mi with rescaled energies 
E™ shifted so that the groundstate energy of H m is zero. 
In the FDM approach, information from different shells 
is combined. This requires that energies from different 
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shells be measured relative to a common groundstate en- 
ergy, which is usually taken to be the absolute ground- 
state energy of the longest chain diagonalized. Hence, 
it is important to keep track of rescaled groundstate en- 
ergies of the H m so that the EJ™ can be related to the 
absolute energies E™ used in the FDM expressions for 
thermodynamic averages (this relation is specific to pre- 
cisely how the sequence H m , m = 1, 2, . . . is defined, so 
we do not specify it here). 

The specific heat C imp (T) = C{T) - C (T) is obtained 
from separate calculations for H and Hq. For H, we first 
calculate E = (H) using (13): 



(H) P = 



N 

E ' 

m—m,Q : l 



,E] n B 



(14) 



I, results in 

kvTxs 



2 = <3> = 



N 

E 

m—rriQ ,S Z ,1 



N 



= J2 hiS)w m B^ 

m—rriQ,l 



(16) 



where £^ S 2 = h(S) = (25 + 1)((2S + l) 2 - 1)/12 has 
been used. For the term xe, we have from (12) 



k B TxE 



(S 



2 ) 

2 . e I 



N 

E 

m—mo,S z ,l,e 



-PET 



w rnS z e - 



(17) 



Since Z m = d N - m Z m and denoting by Z E = d N ~ m the 
partition function of the N — m environment degrees of 
freedom, we can rewrite the above as 



and then substituting this into equation (1) to obtain 



C(T) = k B f3 2 ((H-(H)) 2 ) 

N 

= k B f3 2 ]T w m (Er-E) 2 Br, (15) 

m—rtiQ ,/ 



with a similar calculation for H to obtain Co(T). The 
specific heats are then z-averaged and subtracted to yield 
C imp (T) = C(T)-C (T). Alternatively, the specific heat 
C(T) may be obtained from the z-averaged entropy S(T) 
via C(T) = -TdS/dT, where S is calculated from E 
and Z using S = -dQ/dT = k B InZ + E/T (with sim- 
ilar expressions for Cq(T) and So(T)). We note that in 
cases where explicit numerical derivatives of the thermo- 
dynamic potential with respect to magnetic field or tem- 
perature are required, the NRG supplies a sufficiently 
smooth il(T,B) for this to be possible (see Ref. 12 for 
an early application). We show this within the FDM ap- 
proach for the case of the local magnetic susceptibility 
in Sec. IV C, a quantity that requires a numerical second 
derivative of 0(T, B) with respect to B. 

The susceptibility from (3-4) requires more care, since 
a uniform field acts also on the environment degrees of 
freedom, implying that we require the expectation value 
{{S z + S z , e ) 2 ) in evaluating k B T X (T,B = 0)/(g^ B ) 2 (and 
similarly for k B Txo{T, B — 0)/(g^iB) 2 ), where S z refers 
to total z-component of spin for the system H m and S Zje , 
the total z-component of the N — m environment states e. 
Now {(S z +S z , e ) 2 ) = ((S 2 + S 2 J) since the trace over S z 
(or 5 Zj6 ) of the cross-term 2S z S z ^ e will vanish. Hence, the 
susceptibility will have an additional contribution X e — 
f3{gn B ) 2 (S z e ) due the environment degrees of freedom 
in addition to the usual term xs = P{9^b) 2 {S 2 ) for the 
system H m . Evaluating the latter via (13), indicating 
explicitly the conserved quantum number S z in the trace 
with all other conserved quantum numbers indicated by 
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N 

^ w m Tr 

m-mo,S 2 ,1 

E "mhW^Bi 
m—mo,l 



where f2(S) — J2s = + 1) an d we used the fact 
that TreS* 2 e /Z e = (N — m)/8 since for one environment 
state (Sl e ) ei = Tr ei Sl e ./Z ei = (l/4+l/4)/4 = 1/8 (for 
d = 4). Hence, we have 



k B T X E sr^ f (o\ N ' 



-w m B?, (18) 



m— mo 



and the total susceptibility x{T) 
given by 



XS (T) + X e(T) is 
k B T X (T) _ k B T Xs (T) , k B T XB (T) 



N 

E (h(s) + h(s) 

m—m,Q ,1 



N — m 



w m Br\19) 



with a similar expression for the host susceptibility 
Xo(T). The impurity contribution is then obtained via 
Ximp(T) = x(T) - Xo(T). 

Figure 1 illustrates the problem just discussed for the 
symmetric Anderson model in the strong correlation limit 
(U/A = 12 > 1). Denoting by xs,im P and XE,im P the 
impurity contributions to xs and xej i-e. with respective 
host contributions subtracted, we have Ximp = Xs,imp ~t~ 
XE,imp- We see from Fig. 1 that the contribution from 
the environment degrees of freedom, XE,imp, is significant 
at all temperatures and is required in order to recover the 
exact Bethe ansatz result for Ximp- 
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FIG. 1. (Color online) Impurity contribution to the suscep- 
tibility from Bethe ansatz (xba) and FDM (ximp) vs T/Tk for 
the symmetric Anderson model [U/Ao = 12 and Ao = 0.001D 
with Tk defined in Eq. (20)] . Also shown are the impurity con- 
tributions fc B Txs,im P (T)/(5MB) 2 and kBTxE,in W (T)/(g[i B ) 2 
as defined at the end of Sec. III. The calculations are for 
A = 10 with an energy cut-off e c (A) = 15\/A w 47, with 
z-averaging [n x = 4, z = 1/8, 1/2, 3/8, 3/4]. 



IV. RESULTS 



FIG. 2. (Color online) (a) Impurity specific heat, Ci mp (T), 
and, (b), impurity entropy, Si mp (T), vs reduced tempera- 
ture T/Tk for the symmetric Anderson model with U/Aq = 
12,10,8,6,4,2,1 and A = 0.001D. Broken lines: FDM ap- 
proach. Solid lines: conventional approach. Selected Bethe 
ansatz results are shown as symbols for U /Ao = 12, 8, 4 (cir- 
cles, squares, diamonds, respectively). The Kondo scale is 
defined in Eq. (20). As a guide to the eye, note that the high 
temperature peak in Ci mp shifts downwards with decreasing 
U. NRG and z-averaging parameters as in Fig. 1. 



In this section we compare results for the impurity 
specific heat (Sec. IV A), impurity susceptibilities in re- 
sponse to uniform (Sec. IV B) and local (Sec. IV C) mag- 
netic fields, and the double occupancy (Sec. IV D) of 
the Anderson model, calculated within the FDM ap- 
proach, with corresponding results from the conventional 
approach. For the first two quantities we also show com- 
parisons with Bethe ansatz calculations. Results for the 
Wilson ratio, as a function of Coulomb interaction and 
local level position, within FDM and Bethe ansatz, are 
also presented (Sec. IV B). We show the results for all 
quantities as functions of the reduced temperature T/Tk, 
where the Kondo scale Tk is chosen to be the symmetric 
Anderson model Kondo scale given by 

T K = V^A u 72e- 7rC/ / 8An+7rA °/ 2C/ , (20) 

except for U/Aq < 1 when we set Tk = Ao- The 
Kondo scale in Eq. (20) is related to the T = Bethe 
ansatz susceptibility Xim P (0) via Xim P (0) = (gAi B ) 2 /4Tk 
in the limit U 3> Ao (see Ref. 7). We continue to 
use Tk in comparing results from different methods, al- 
though we note that for the asymmetric Anderson model, 
the physical low energy Kondo scale, Tl, will increas- 
ingly deviate from Tk with increasing level asymmetry 
S = 2sd + U. For example, second order poor Man's scal- 
ing for the Anderson model yields a low energy scale' 55 
T L = y/UA^j2e ned( - ed+u ^ 2AaU . 



A. Specific heat 



Figure 2 shows the impurity specific heat (Ci mp ) and 
impurity entropy (Si mp ) for the symmetric Anderson 
model versus temperature T/Tk for increasing Coulomb 
interaction U/Aq. One sees from Fig. 2 that there is 
excellent agreement between the results obtained within 
the FDM approach, within the conventional approach 
and within the exact Bethe ansatz calculations. This 
agreement is found for both the strongly correlated limit 
U/Aq 1 where there are two peaks in the specific 
heat, a low temperature Kondo induced peak and a high 
temperature peak due to the resonant level, and for the 
weakly correlated limit U/Aq < 1, where there is only 
a single resonant level peak in the specific heat. The 
correct high temperature entropy In 4 is obtained in all 
cases. 

The temperature dependence of the impurity specific 
heat and entropy, for the asymmetric Anderson model 
is shown in Fig. 3 for local level positions ranging from 
£rf/A = —6 to +5 in units of A . For simplicity we 
continue to show the results as a function of T/Tk, with 
Tk the symmetric Kondo scale (20), although, the true 
Kondo scale will deviate from this for Ed > —U/2. The 
FDM results agree also here very well with the conven- 
tional approach and the Bethe ansatz calculations. 



FIG. 3. (Color online) (a) Impurity specific heat, Ci mp (T), 
and, (b), impurity entropy Si mp (T), vs reduced temperature 
T/Tk for the asymmetric Anderson model with U/Ao = 12, 
Ao = 0.001D and several values of Ed/ Ao = —5, —3, . . . , +5. 
Broken lines: FDM approach. Solid lines: conventional ap- 
proach. Bethe ansatz results are shown as symbols for Ed/A. 
For simplicity we used the symmetric Tk of Eq. (20) for all 
Ed values. NRG and z-averaging parameters as in Fig. 1. 



FIG. 4. (Color online) Impurity susceptibility, Ximp(T), 
vs T/Tk for the symmetric Anderson model with U /Ao = 
12,10,8,6,4,2,1,0.01 and A = 0.001D with T K defined in 
Eq. (20) for U/A > 1 and T K = A for the U/A = 0.01 
case. Broken lines: FDM approach. Solid lines: conven- 
tional approach. Symbols: Bethe ansatz (for selected values 
of U/Ao = 12,8,6,4,2). As a guide to the eye, note that 
Ximp is increasingly enhanced with increasing U. NRG and 
z-averaging parameters as in Fig. 1. 
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0.2637 


1.852 


1.877 


2 




0.3085 




1.578 


1 




0.2214 




1.317 


0.01 




0.1599 




1.003 



TABLE I. Zero temperature susceptibilities 
kvTKXim-p/igw) 2 and Wilson ratios R a = 
lim T ^o47r 2 xf mp (T')/3Cf m p(r)/r for the symmetric An- 
derson model at several values of U/Ao using the Bethe 
ansatz/NRG FDM approach (a = B A/NRG). Note that T K 
is defined by Eq. 20 for U/Ao > 1 and is set to Ao otherwise. 



B. Susceptibility and Wilson ratio 

Fig. 4 compares the susceptibilities of the symmet- 
ric Anderson model calculated from FDM, conven- 
tional and Bethe ansatz approaches for several val- 
ues of U/Aq, again indicating good agreement over 
the whole temperature range between these three ap- 
proaches. Table I lists the zero temperature impurity 
susceptibilities (kBTxXinip/igHB) 2 ) and Wilson ratios 
(R = lim T ^o 47r 2 v im p(T)/3C im p(T)/T), for the symmet- 
ric Anderson model as calculated within FDM and for 
a range of Coulomb interactions from strong U / A 1 
to weak [U/Aq <C 1). In these two limits, the Wilson 
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FIG. 5. (Color online) Impurity susceptibility, Xim P (T), vs 
T/Tk for the asymmetric Anderson model with U/Ao = 12, 
Ao = 0.001-D and several values of £d/Ao with Tk defined in 
Eq. (20). Broken lines: FDM approach. Solid lines: conven- 
tional approach. Symbols: Bethe ansatz (for selected values 
of Ed/ Ao — — 5, — 3, — 1, 0, +1, +3). As a guide to the eye, 
note that the susceptibility curves shift to higher tempera- 
tures with increasing Ed- NRG and z-averaging parameters 
as in Fig. 1. 

ratio for the symmetric Anderson model approaches the 
well known values of 2, and 1, respectively within FDM 
(a = NRG) and Bethe ansatz. Comparison with Bethe 
ansatz results at selected values of U/Aq indicate an er- 
ror in the susceptibility of around 2% with a similar error 
in the Wilson ratio. 



7 



Ed/ t-±Q 


kBT K Ximp 


/(SMs) 2 


R a 




a = BA 


a = NRG 


a = BA a 


= NRG 


-5 


0.219482 


0.2245 


1.999 


2.025 


-4 


- 


0.1515 


- 


2.023 


-3 


0.077356 


0.0785 


1.990 


2.00 


-2 


- 


0.0315 


- 


1.97 


-1 


0.010337 


0.0103 


1.795 


1.78 





0.003303 


0.0033 


1.512 


1.50 


1 
1 


U.UU1ZOU 


n nni ^ 

U.UUlo 


i.oio 




2 




0.00059 




1.18 


3 


0.000325 


0.00033 


1.086 


1.12 


4 




0.00021 




1.09 


5 




0.00014 




1.06 



TABLE II. Zero temperature susceptibilities fcBTx?m P /(3MB) 2 
and Wilson ratios R a = lim T ^ 47r 2 xf mp (T)/3C' i a mp (T)/T for 
the asymmetric Anderson model at U /Ao = 12 and several 
local level positions Ed/ Ao using the Bethe ansatz/FDM NRG 
approach (a = BA/NRG). 



Fig. 5 shows results within FDM and conventional 
approaches for the asymmetric Anderson model (ed > 
— {7/2) and for several local level positions ranging from 
the Kondo (— £d/Ao ~> 1) to the mixed valence \ed/Ao\ < 
1 and empty orbital regimes e,j/Ao > 1. Bethe ansatz 
results are also shown for selected local level positions, 
and we see again very good agreement between all three 
methods over the whole temperature range. Correspond- 
ing zero temperature susceptibilities and Wilson ratios 
are listed in Table II. Note that the Wilson ratio ap- 
proaches the value for a non-interacting system only in 
the empty orbital limit (e^ 3> Ao), being approximately 
1.5 ± 0.25 in the mixed valence regime (|£d/A | < 1). 
The Wilson ratio from NRG and Bethe ansatz deviate 
by less than 3% in all regimes. 



C. Local susceptibility 

It is also interesting to consider the susceptibility, xioc, 
in response to a local magnetic field acting only at the 
impurity site and to compare this with the susceptibil- 
ity, Ximp, discussed above, in which the magnetic field 
acts on both the impurity and conduction electron spins. 
The former is relevant, for example, in nuclear magnetic 
resonance and neutron scattering experiments, while the 
latter can be measured in bulk samples with and without 
magnetic impurities. 

A local magnetic field term, — g/iB-B^ ^, in the An- 
derson model, with S Xi d = (n^t — n d^)/2, is not a con- 
served quantity, i.e S z< d is not conserved , and xioc(T) 
cannot be expressed as a fluctuation as in Eq. (3-4), 
which would obviate the need to explicitly evaluate a 
numerical second derivative with respect to B of the 
thermodynamic potential. Such a derivative, however, 
poses no actual problem within NRG, so we proceed by 
explicitly diagonalizing the Anderson model in a local 
field, using only U(l) symmetries for charge and spin 




10" 3 10" 2 



FIG. 6. ( Color online) Comparison of the local, Xioc(r), and 
impurity, Ximp, susceptibilities vs T/Tk for the symmetric An- 
derson model with U/A = 12,8,4,2 and A = 0.001D with 
7k defined in Eq. (20). Broken lines: FDM approach. Sym- 
bols: Bethe ansatz (for selected values of U/Ao = 12,8,4,2). 
NRG and z-averaging parameters as in Fig. 1. 



(for the symmetric Anderson model in a magnetic field, 
an SU(2) pseudo-spin symmetry may be exploited, by 
using the mapping of this model in a local magnetic 
field B onto the SU(2) invariant negative-i7 Anderson 
model in zero magnetic field at finite level asymmetry 
2ed + U = B 3e ' 37 ). The evaluation of xioc then proceeds 
via xioc(T,S = 0) = -d 2 n Xoc {T,B)/dB 2 \ B=0 where 
fl ioc (T,B) = Q(T,B) - n (T) and Sl(T,B) and n (T) 
are the thermodynamic potentials of the total system in 
a local magnetic field B and the host system, respectively. 

Results for xioc obtained in this way are shown in Fig. 6 
at several values of U /Ao as a function of T/Tk- A com- 
parison of xioc to Ximp obtained from the Bethe ansatz, 
allows us to conclude that these two susceptibilities are 
close to identical at all temperatures, i.e. Ximp(7 1 ) = 
Xioc{T) and for all interaction strengths U/A . This 
is not always the case. A prominent example is the 
anisotropic Kondo model 38 , where Ximp = «Xioc, with 
the dissipation strength < a < 1 being determined by 
the anisotropy of the exchange interaction 38 ' 39 . 

Figure 7 compares local and impurity susceptibilities 
for the asymmetric Anderson model in the strong corre- 
lation limit (U/Aq = 12) for several local level positions, 
ranging from the Kondo (e^/Ao = —5, —4, —3, —2) to the 
mixed valence (e^/Ao = — 1,0, +1) and into the empty 
orbital regime (e^/Ao = +2, ...,+5). We see that, as 
for the symmetric Anderson model, local and impurity 
susceptibilities are almost identical at all temperatures 
and for all local level positions, i.e. Ximp(T) — xioc(T) 
for the parameter values used. 

The result Ximp(7 1 ) = Xioc(7 1 ), which we verified 
here, follows from the Clogston- Anderson compensation 
theorem 40 (see Ref. 7). Consider the impurity contribu- 



FIG. 7. (Color online) Comparison of the local, Xioc, and 
impurity, Xim P (T), susceptibilities vs T/Tk for the asymmet- 
ric Anderson model with U/A = 12, A = 0.001D with T K 
defined in Eq. (20) and for several values of the local level 
position: £d/Ao = — 5, — 4, . . . , +5. Broken lines: xioc from 
the FDM approach. Symbols: Xim P from Bethe ansatz for se- 
lected values of Ed- The NRG calculations are for A = 4 with 
an energy cut-off e c (A = 4) = 40 with z-averaging [n z = 2 
with z = and z = 0.5]. 



FIG. 8. (Color online) (a) Double occupancy D occ = 
{ndt n di) as a function of temperature T/Tk for the sym- 
metric Anderson model and decreasing values of U/Aq = 
12,10,8,6,4,2,1,0.01 within FDM (solid lines) and conven- 
tional approaches (symbols), (b) Double occupancy D occ = 
{nd^ndi) as a function of temperature T/Tk for the asym- 
metric Anderson model and increasing values of Ed/ Ao = 
-5, -3,-1,0, +1, +3, +5 for U/A = 12 within FDM (solid 
lines) and conventional approaches (symbols). Tk is defined 
in Eq. (20) for U /Ao > 1 and is set to Ao for the case 
U/Aq = 0.01. NRG and z-averaging parameters as in Fig. 1. 



tion to the magnetization Mi mp (B) in a uniform field. 
This is given by M imp /(gfj, B ) = (S z ,d + S z<c ) - (S ZtC } 
where S z .d, S ZtC are the impurity and conduction electron 
z— components of spin. Using equations of motion, one 
easily shows ' , that the additional impurity magnetiza- 
tion 5Mi mp / (<?/Ub) = (S z ,c) — (S z .c)o from the conduction 
electrons induced by the presence of the impurity is given 
by 



SM h 



(#Mb) 



2^5Z / '^'/Uilm 



aG da {u,B) 



dA(u) 



du 



V) 



where f(u>) is the Fermi function, Gda(uJ, B) is the spin 
a local level Green function of the Anderson model 
and A(lj) = J2k \Vk\ 2 /(u ~ £fe + iS) is the hybridiza- 
tion function. For a flat band, dA(uj)/du w Aq/D 
in the wide-band limit. Hence, 5Mi mp /(gfj,B) is of or- 
der {M loc (B)/g[x B )A /D where M loc (B)/(g^ B ) = (n df - 
Tidi) /2 = (S z ^d) is the local magnetization, which is linear 
in B for B — > 0. From this we deduce that Mi mp /(g/j, B ) ss 
(S z d) = Mi oc /(gfi B ) to within corrections of order 
(M loc (B)/ gf i B )A Q /D, i.e. (T) = Xioc(T) to within 
corrections of order (x\ oc (T)/(giJ, B ) 2 )A a / D < X\oc(T). 
Away from the wide-band limit, or for strong energy de- 
pendence of A(cj), the above susceptibilities will differ by 
the correction term given by the field derivative of 5M lmp 
in Eq. (21). 



D. Double occupancy 

Our conclusions concerning the accuracy of specific 
heat and the susceptibility calculations within the FDM 
approach, hold also for other thermodynamic properties, 
e.g. for the occupation number or the double occupancy. 
Fig. 8(a) shows a comparison between the FDM and con- 
ventional approaches for the temperature dependence of 
the double occupancy D occ = (ndfUdi) of the sym- 
metric model for different strengths of correlation U/A , 
and Fig. 8(b) shows the same for the asymmetric An- 
derson model for U/Ao = 12 and for several local level 
positions. The results of the two approaches agree at all 
temperatures, local level positions and Coulomb inter- 
actions. Notice that D acquires its mean-field value of 
1/4 for the symmetric model in the limit U/Aq — > and 
is strongly suppressed with increasing Coulomb interac- 
tion away from this limit [see Fig. 8(a)]. Similarly for 
the asymmetric model, increasing Ed/Ao away from the 
correlated Kondo regime decreases the double occupancy 
significantly [see Fig. 8(b)]. 



V. SUMMARY 

In this paper we focused on the calculation of the im- 
purity specific heat and the impurity susceptibility of 
the Anderson model within the FDM approach 1 , finding 
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that this method gives reliable results for these quan- 
tities, as shown by a comparison to both exact Bethe 
ansatz calculations 22-28 and to NRG calculations within 
the conventional approach 29 . Some care is needed in im- 
plementing the FDM approach for the susceptibility Ximp 
in a uniform field, i.e. when the applied magnetic field 
acts on both the impurity and conduction electron spins. 
In this case, an additional contribution from the envi- 
ronment degrees of freedom needs to be included. We 
also showed that the susceptibility in response to a local 
magnetic field on the impurity, X\oa could also be ob- 
tained within FDM and a comparison of this susceptibil- 
ity with Ximp (from Bethe ansatz) , showed that they are 
close to identical at all temperatures, and in all parame- 
ter regimes for A <C D, thereby verifying the Clogston- 
Anderson compensation theorem. An arbitrary temper- 
ature grid can be used for thermodynamics in both the 
conventional and the FDM approaches, however, the for- 
mer requires a specific best shell to be selected depending 
on T and z, whereas the FDM approach avoids this step 



by incorporating all excitations from all shells in a single 
density matrix. 

We also showed, that quantities such as the double 
occupancy can also be accurately calculated within the 
FDM approach. The double occupancy can be probed in 
experiments on cold atom realizations of Hubbard mod- 
els in optical lattices. 41,42 Flexible techniques, such as the 
FDM approach, for calculating them within a dynamical 
mean field theory 4 treatment of the underlying effec- 
tive quantum impurity models could be useful in future 
investigations of such systems 46 . 

ACKNOWLEDGMENTS 

We thank Jan von Delft, Markus Hani and Ralf Bulla 
for useful discussions and acknowledge supercomputer 
support by the John von Neumann institute for Comput- 
ing (Jiilich). Support from the DFG under grant number 
WE4819/1-1 is also acknowledged (AW). 



1 A. Weichselbaum and J. von Delft, Phys. Rev. Lett. 99, 
076402 (2007). 

2 F. B. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 
(2005). 

3 K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975). 

4 H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson, 
Phys. Rev. B 21, 1003 (1980). 

5 H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson, 
Phys. Rev. B 21, 1044 (1980). 

6 R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 
80, 395 (2008). 

7 A. C. Hewson, The Kondo Problem to Heavy Fermions 
(Cambridge University Press, Cambridge, 1997). 

8 L. N. Oliveira and J. W. Wilkins, Phys. Rev. Lett. 47, 
1553 (1981). 

9 H. O. Frota and L. N. Oliveira, Phys. Rev. B 33, 7871 
(1986). 

10 O. Sakai, Y. Shimizu, and T. Kasuya, Journal of the Phys- 
ical Society of Japan 58, 3666 (1989). 

11 T. A. Costi and A. C. Hewson, Philosophical Magazine 
Part B 65, 1165 (1992). 

12 T. A. Costi, A. C. Hewson, and V. Zlatic, J. Phys.: Con- 
dens. Matter 6, 2519 (1994). 

13 R. Bulla, A. C. Hewson, and T. Pruschke, Journal of 
Physics: Condensed Matter 10, 8365 (1998). 

14 W. Hofstetter, Phys. Rev. Lett. 85, 1508 (2000). 

15 T. A. Costi, Phys. Rev. B 55, 3003 (1997). 

16 F. B. Anders, Phys. Rev. Lett. 101, 066804 (2008). 

17 R. Peters, T. Pruschke, and F. B. Anders, Phys. Rev. B 
74, 245114 (2006). 

18 A. I. Toth, C. P. Moca, O. Legeza, and G. Zarand, Phys. 
Rev. B 78, 245109 (2008). 

19 T. A. Costi, L. Bergqvist, A. Weichselbaum, J. von Delft, 
T. Micklitz, A. Rosch, P. Mavropoulos, P. H. Dederichs, 
F. Mallet, L. Saminadayar, and C. Bauerle, Phys. Rev. 
Lett. 102, 056802 (2009). 

20 T. A. Costi and V. Zlatic, Phys. Rev. B 81, 235127 (2010). 

21 P. W. Anderson, Physical Review 124, 41 (1961). 



22 N. Kawakami and A. Okiji, Physics Letters A 86, 483 
(1981). 

23 N. Kawakami and A. Okiji, Journal of the Physical Society 
of Japan 51, 2043 (1982). 

24 A. Okiji and N. Kawakami, Phys. Rev. Lett. 50, 1157 
(1983). 

25 P. B. Wiegmann and A. M. Tsvelick, Journal of Physics 
C: Solid State Physics 16, 2281 (1983). 

26 A. M. Tsvelick and P. B. Wiegmann, Journal of Physics 
C: Solid State Physics 16, 2321 (1983). 

27 V. M. Filyov, A. M. Tsvelick, and P. B. Wiegmann, 
Physics Letters A 89, 157 (1982). 

28 A. M. Tsvelick and P. B. Wiegmann, Physics Letters A 89, 
368 (1982). 

29 V. L. Campo and L. N. Oliveira, Phys. Rev. B 72, 104432 
(2005). 

30 L. Merker and T. A. Costi, Phys. Rev. B 86, 075150 (2012). 

31 W. C. Oliveira and L. N. Oliveira, Phys. Rev. B 49, 11986 
(1994). 

32 R. Zitko and T. Pruschke, Phys. Rev. B 79, 085106 (2009). 

33 S. C. Costa, C. A. Paula, V. L. Lfbero, and L. N. Oliveira, 
Phys. Rev. B 55, 30 (1997). 

34 A. Weichselbaum, Phys. Rev. B 84, 125130 (2011). 

35 F. D. M. Haldane, Phys. Rev. Lett. 40, 416 (1978). 

36 G. Iche and A. Zawadowski, Solid State Communications 
10, 1001 (1972). 

37 A. C. Hewson, J. Bauer, and W. Roller, Phys. Rev. B 73, 
045117 (2006). 

38 P. B. Vigman and A. M. Finkelstein, Sov. Phys. JETP 48, 
102 (1978). 

39 F. Guinea, V. Hakim, and A. Muramatsu, Phys. Rev. B 
32, 4410 (1985). 

40 A. M. Clogston and P. W. Anderson, Bull. Am. Phys. Soc. 
6, 124 (1961). 

41 R. Jordens, N. Strohmaier, K. Giinter, H. Moritz, and 
T. Esslinger, Nature (London) 455, 204 (2008). 

42 U. Schneider, L. Hackermuller, S. Will, T. Best, I. Bloch, 
T. A. Costi, R. W. Helmes, D. Rasch, and A. Rosch, 



10 



Science 322, 1520 (2008). 

G. Kotliar and D. Vollhardt, Physics Today 57, 53 (2004). 
A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, 
Rev. Mod. Phys. 68, 13 (1996). 



D. Vollhardt, Annalen der Physik 524, 1 (2012). 

B. Tang, T. Paiva, E. Khatami, and M. Rigol, ArXiv e- 

prints (2012), arXiv:1206.0006 [cond-mat.str-el]. 



